### Beginning ####
rm(list = ls())
# analysis study 1
pacman::p_load(multiwayvcov,
               margins, 
               haven, 
               stargazer,
               dplyr, 
               lmerTest, 
               lmtest, 
               fastDummies, 
               purrr, 
               tidyr, 
               stringr,
               huxtable, 
               parallel, 
               modelsummary, 
               readr
               )
# devtools::install_github("ChandlerLutz/starpolishr")
library(starpolishr)

psych_long <- read_dta("./data_clean/psych_clean_long.dta") %>% 
  rename(StoryTrue = true,
         CRT = zCRT_all,
         CRTUnscaled = CRT_all,
         # AOT = AOT_all,
         AOTUnscaled =  AOT_allUnscaled,
         # AOT_shortUnscaled,
         # AOT_longUnscaled,
         RussianOrientation = russianOrientation,
         ProRussiaPoliOrientation = zantiEuropeOrientation,
         RussianOrientationUnscaled = XrussianOrientation,
         `78.5` = antiEuropeOrientation,
         Female = RSPSEX,
         Age = RSPAGE,
         Education = RSPEDUC,
         Income = income) %>% 
  mutate(Belief = response,
         AOT = (AOTUnscaled - mean(AOTUnscaled)) / sd(AOTUnscaled),
         AOTUnscaled = (AOTUnscaled - 1)/4)

psych_short <- read_dta("./data_clean/psych_clean_short.dta") %>% 
  rename(
    # StoryTrue = true,
    CRT = zCRT_all, #z-scaled
    CRTUnscaled = CRT_all,
    # AOT = AOT_all,
    AOTUnscaled =  AOT_allUnscaled,
    RussianOrientation = russianOrientation, #also z-scaled, variable with x contains unscaled
    ProRussiaPoliOrientation = zantiEuropeOrientation, #z-scaled
    RussianOrientationUnscaled = XrussianOrientation,
    ProRussiaPoliOrientationUnscaled = antiEuropeOrientation,
    Female = RSPSEX,
    Age = RSPAGE,
    Education = RSPEDUC,
    Income = income) %>% 
  mutate(AOT = (AOTUnscaled - mean(AOTUnscaled)) / sd(AOTUnscaled),
         AOTUnscaled = (AOTUnscaled - 1)/4)

# +++++++++++++++++++++++++++++
# Summary Statistics Table ####
# +++++++++++++++++++++++++++++

sum_funs <- list(mean = function(x) mean(x, na.rm = TRUE),
                 median = function(x) median(x, na.rm = TRUE), 
                 min = function(x)min(x, na.rm = TRUE),
                 max = function(x)max(x, na.rm = TRUE),
                 se = function(x)sd(x, na.rm = TRUE)/sqrt(n()), 
                 sd = function(x)sd(x, na.rm = TRUE))

# scaled variables

df.sum <- cbind(
  summarise_all(
    select(psych_short, CRT, AOT, RussianOrientation, ProRussiaPoliOrientation, Female, Age, Education, Income),
    .funs = sum_funs),
  summarise_all(
    select(psych_long, Belief, StoryTrue), .funs = sum_funs
  )
) 

df.stats.tidy <- 
  df.sum %>% gather(stat, val) %>%
  separate(stat, into = c("var", "stat"), sep = "_") %>%
  spread(stat, val) %>% 
  mutate(var = factor(var, levels = c("Belief", "StoryTrue",
                                      "CRT", "AOT", 
                                      "RussianOrientation",
                                      "ProRussiaPoliOrientation",
                                      "Education", "Income", "Age", "Female"))) %>% 
  select(var, min, mean, median, max, se, sd) %>% # reorder columns
  rename(Variable = var,
         Min. = min,
         Mean = mean,
         Median = median,
         Max. = max, 
         `Std. Error` = se,
         `Std. Dev.` = sd); # rename columns

df.stats.tidy <- df.stats.tidy[order(df.stats.tidy$Variable), ] # reorder manually by factor levels for Variable

t <- xtable::xtable(df.stats.tidy, type = "latex")
print(t, file = "./tables/sumstats_s1.tex", include.rownames = FALSE)

# unscaled variables

df.sum_us <- cbind(
  summarise_all(
    select(psych_short, CRTUnscaled, AOTUnscaled, RussianOrientationUnscaled, ProRussiaPoliOrientationUnscaled, 
           Female, Age, Education, Income),
    .funs = sum_funs),
  summarise_all(
    select(psych_long, Belief, StoryTrue), .funs = sum_funs
  )
)

df.stats.tidy_us <- 
  df.sum_us %>% gather(stat, val) %>%
  separate(stat, into = c("var", "stat"), sep = "_") %>%
  spread(stat, val) %>% 
  mutate(var = factor(var, levels = c("Belief", "StoryTrue",
                                      "CRTUnscaled", "AOTUnscaled", 
                                      "RussianOrientationUnscaled",
                                      "ProRussiaPoliOrientationUnscaled",
                                      "Education", "Income", "Age", "Female"))) %>% 
  select(var, min, mean, median, max, se, sd) %>% # reorder columns
  rename(Variable = var,
         Min. = min,
         Mean = mean,
         Median = median,
         Max. = max, 
         `Std. Error` = se,
         `Std. Dev.` = sd); # rename columns

df.stats.tidy_us <- df.stats.tidy_us[order(df.stats.tidy_us$Variable), ] # reorder manually by factor levels for Variable

t_us <- xtable::xtable(df.stats.tidy_us, type = "latex")
print(t_us, file = "./tables/sumstats_s1_unscaled.tex", include.rownames = FALSE)

# ++++++++++++++++++++++++++++
# Zero-Order Correlations ####
# ++++++++++++++++++++++++++++

# full correlations CRT, AOT, and each of the 16 forms of statement

# map function onto a split - need cor between outcome and each one of the 16 cells
# show first few columns, no need for whole table

psych_long <-
  psych_long %>% mutate(
    strat_by_theme = case_when(
      theme == "economic" & strategy == 1 ~ "economic|antiwest",
      theme == "economic" & strategy == 2 ~ "economic|proRussia",
      theme == "economic" & strategy == 3 ~ "economic|undermineUkraine",
      theme == "economic" & strategy == 4 ~ "economic|truestory",
      theme == "political" & strategy == 1 ~ "political|antiwest",
      theme == "political" & strategy == 2 ~ "political|proRussia",
      theme == "political" & strategy == 3 ~ "political|undermineUkraine",
      theme == "political" & strategy == 4 ~ "political|truestory",
      theme == "historical" & strategy == 1 ~ "historical|antiwest",
      theme == "historical" & strategy == 2 ~ "historical|proRussia",
      theme == "historical" & strategy == 3 ~ "historical|undermineUkraine",
      theme == "historical" & strategy == 4 ~ "historical|truestory",
      theme == "military" & strategy == 1 ~ "military|antiwest",
      theme == "military" & strategy == 2 ~ "military|proRussia",
      theme == "military" & strategy == 3 ~ "military|undermineUkraine",
      theme == "military" & strategy == 4 ~ "military|truestory",
    )
  )

# Full Sample
cortabs <- 
  psych_long %>% 
  split(.$strat_by_theme) %>% 
  map(~ cor(select(., CRT, AOT, response))) %>% 
  data.frame() %>%
  mutate(X = rownames(.)) %>% 
  pivot_longer(-X) %>% 
  separate(name, c("Topic", "Strategy", "Y"), "[.]") %>% 
  filter(value != 1) %>%
  mutate(value = round(value, 3),
         Strategy = factor(Strategy, levels = c("antiwest", "proRussia", "undermineUkraine", "truestory"))) %>% 
  arrange(X, Topic, Strategy) %>%  
  pivot_wider(names_from = Y, values_from = value) %>% 
  split(.$X) %>% 
  lapply(., "[", c(-1))

cortabs$AOT <- cortabs$AOT[, c(-3, -5)]
cortabs$CRT <- cortabs$CRT[, c(-3, -5)]
cortabs$response <- cortabs$response[, -4]

stargazer(cortabs$response %>% 
            mutate(Strategy = as.character(Strategy)), 
          type = "latex", 
          title = "Study 1 Correlations with Belief in Story",
          summary = FALSE, 
          header = FALSE,
          out = "tables/cors1_s1.tex")

# Pro Europe Subset
cortabs <- 
  psych_long %>% 
  filter(proEurope == 1) %>% 
  split(.$strat_by_theme) %>% 
  map(~ cor(select(., CRT, AOT, response))) %>% 
  data.frame() %>%
  mutate(X = rownames(.)) %>% 
  pivot_longer(-X) %>% 
  separate(name, c("Topic", "Strategy", "Y"), "[.]") %>% 
  filter(value != 1) %>%
  mutate(value = round(value, 3),
         Strategy = factor(Strategy, levels = c("antiwest", "proRussia", "undermineUkraine", "truestory"))) %>% 
  arrange(X, Topic, Strategy) %>%  
  pivot_wider(names_from = Y, values_from = value) %>% 
  split(.$X) %>% 
  lapply(., "[", c(-1))

cortabs$AOT <- cortabs$AOT[, c(-3, -5)]
cortabs$CRT <- cortabs$CRT[, c(-3, -5)]
cortabs$response <- cortabs$response[, -4]

stargazer(cortabs$response %>% 
            mutate(Strategy = as.character(Strategy)), 
          type = "latex", 
          title = "Study 1 Correlations with Belief in Story (Pro-Europe)",
          summary = FALSE, 
          header = FALSE,
          out = "tables/cors2_s1.tex")

# Anti Europe Subset
cortabs <- 
  psych_long %>% 
  filter(antiEurope == 1) %>% 
  split(.$strat_by_theme) %>% 
  map(~ cor(select(., CRT, AOT, response))) %>% 
  data.frame() %>%
  mutate(X = rownames(.)) %>% 
  pivot_longer(-X) %>% 
  separate(name, c("Topic", "Strategy", "Y"), "[.]") %>% 
  filter(value != 1) %>%
  mutate(value = round(value, 3),
         Strategy = factor(Strategy, levels = c("antiwest", "proRussia", "undermineUkraine", "truestory"))) %>% 
  arrange(X, Topic, Strategy) %>%  
  pivot_wider(names_from = Y, values_from = value) %>% 
  split(.$X) %>% 
  lapply(., "[", c(-1))

cortabs$AOT <- cortabs$AOT[, c(-3, -5)]
cortabs$CRT <- cortabs$CRT[, c(-3, -5)]
cortabs$response <- cortabs$response[, -4]

stargazer(cortabs$response %>% 
            mutate(Strategy = as.character(Strategy)), 
          type = "latex", 
          title = "Study 1 Correlations with Belief in Story (Anti-Europe)",
          summary = FALSE, 
          header = FALSE,
          out = "tables/cors3_s1.tex")

# Far Right Subset
cortabs <- 
  psych_long %>% 
  filter(farRight == 1) %>% 
  split(.$strat_by_theme) %>% 
  map(~ cor(select(., CRT, AOT, response))) %>% 
  data.frame() %>%
  mutate(X = rownames(.)) %>% 
  pivot_longer(-X) %>% 
  separate(name, c("Topic", "Strategy", "Y"), "[.]") %>% 
  filter(value != 1) %>%
  mutate(value = round(value, 3),
         Strategy = factor(Strategy, levels = c("antiwest", "proRussia", "undermineUkraine", "truestory"))) %>% 
  arrange(X, Topic, Strategy) %>%  
  pivot_wider(names_from = Y, values_from = value) %>% 
  split(.$X) %>% 
  lapply(., "[", c(-1))

cortabs$AOT <- cortabs$AOT[, c(-3, -5)]
cortabs$CRT <- cortabs$CRT[, c(-3, -5)]
cortabs$response <- cortabs$response[, -4]

stargazer(cortabs$response %>% 
            mutate(Strategy = as.character(Strategy)), 
          type = "latex", 
          title = "Study 1 Correlations with Belief in Story (Far-Right)",
          summary = FALSE, 
          header = FALSE,
          out = "tables/cors4_s1.tex")

# correlations CRT, AOT, 
#support for Ukraine vs Russia,
# membership in Pro-Europe vs Pro-Russia, Far Right vs Other, 
# having voted, income, education

cordat <- 
  psych_short %>% 
  dplyr::select(
    CRT, 
    AOT, 
    RussianOrientation,
    ProRussiaPoliOrientation,
    proEurope,
    farRight,
    VOTPRELE,
    Income,
    Education)

cordat %>% 
  cor(use = "complete.obs") %>% 
  stargazer(font.size = "scriptsize", 
            title = "Study 1: Zero-Order Correlations with Covariates",
            out = "tables/cors5_s1.tex",
            type = "latex")

# +++++++++++++++++++++++++++
# Linear Models - scaled ####
# +++++++++++++++++++++++++++

# part 1
m1_s1 <- list()
m1_s1[[1]] <- lm(response ~ CRT*StoryTrue, psych_long)
m1_s1[[2]] <- lm(response ~ AOT*StoryTrue, psych_long)
m1_s1[[3]] <- lm(response ~ cogStyleAgg*StoryTrue, psych_long)

vcov1 <- list()
vcov1 <- lapply(m1_s1, cluster.vcov, ~ Respondent_Serial + question)

# part 2
m2_s1 <- list()
m2_s1[[1]] <- lm(response ~ CRT*StoryTrue*RussianOrientation, psych_long)
m2_s1[[2]] <- lm(response ~ AOT*StoryTrue*RussianOrientation, psych_long)
m2_s1[[3]] <- lm(response ~ CRT*StoryTrue*ProRussiaPoliOrientation, psych_long)
m2_s1[[4]] <- lm(response ~ AOT*StoryTrue*ProRussiaPoliOrientation, psych_long)

vcov2 <- list()
vcov2 <- lapply(m2_s1, cluster.vcov, ~ Respondent_Serial + question)

# part 3
m3_s1 <- list()
m3_s1[[1]] <- lm(response ~ CRT*StoryTrue*RussianOrientation + Age + Female + factor(Education) + Income, psych_long)
m3_s1[[2]] <- lm(response ~ AOT*StoryTrue*RussianOrientation + Age + Female + factor(Education) + Income, psych_long)
m3_s1[[3]] <- lm(response ~ CRT*StoryTrue*ProRussiaPoliOrientation + Age + Female + factor(Education) + Income, psych_long)
m3_s1[[4]] <- lm(response ~ AOT*StoryTrue*ProRussiaPoliOrientation + Age + Female + factor(Education) + Income, psych_long)

vcov3 <- list()
vcov3 <- lapply(m3_s1, cluster.vcov, ~ Respondent_Serial + question)

# part 4 (with FEs for theme)
m4_s1 <- list()
m4_s1[[1]] <- lm(response ~ CRT*StoryTrue*RussianOrientation + Age + Female + factor(Education) + Income + factor(theme), 
                 psych_long)
m4_s1[[2]] <- lm(response ~ AOT*StoryTrue*RussianOrientation + Age + Female + factor(Education) + Income + factor(theme), 
                 psych_long)
m4_s1[[3]] <- lm(response ~ CRT*StoryTrue*ProRussiaPoliOrientation + Age + Female + factor(Education) + Income + 
                   factor(theme), 
                 psych_long)
m4_s1[[4]] <- lm(response ~ AOT*StoryTrue*ProRussiaPoliOrientation + Age + Female + factor(Education) + Income + 
                   factor(theme), 
                 psych_long)

vcov4 <- list()
vcov4 <- lapply(m4_s1, cluster.vcov, ~ Respondent_Serial + question)

# get SEs
cl.rob.se.1 <- lapply(vcov1, function(x) sqrt(diag(x)))
cl.rob.se.2 <- lapply(vcov2, function(x) sqrt(diag(x)))
cl.rob.se.3 <- lapply(vcov3, function(x) sqrt(diag(x)))

# get F Stats
cl.wald1 <- list()
for(i in 1:length(m1_s1)) {
  cl.wald1[[i]]  <- waldtest(m1_s1[[i]], vcov = vcov1[[i]])
}

cl.wald2 <- list()
for(i in 1:length(m2_s1)) {
  cl.wald2[[i]]  <- waldtest(m2_s1[[i]], vcov = vcov2[[i]])
}

cl.wald3 <- list()
for(i in 1:length(m3_s1)) {
  cl.wald3[[i]]  <- waldtest(m3_s1[[i]], vcov = vcov3[[i]])
}

# get p-values
p.vals1 <- list()
for(i in 1:length(m1_s1)) {
  p.vals1[[i]] <- coeftest(m1_s1[[i]], vcov = vcov1[[i]])[,4]
}

p.vals2 <- list()
for(i in 1:length(m2_s1)) {
  p.vals2[[i]] <- coeftest(m2_s1[[i]], vcov = vcov2[[i]])[,4]
}

p.vals3 <- list()
for(i in 1:length(m3_s1)) {
  p.vals3[[i]] <- coeftest(m3_s1[[i]], vcov = vcov3[[i]])[,4]
}

# saving main models
main_models_study1 <- list(m1_s1, m2_s1, m3_s1)
save(main_models_study1, file =  "./data_clean/main_clustered_models_study1.Rdata")

load("./data_clean/main_clustered_models_study1.Rdata")

# exporting model results

mod.names <- vector("character", sum(length(m1_s1), length(m2_s1), length(m3_s1)))
for(i in 1:sum(length(m1_s1), length(m2_s1), length(m3_s1))) {
  mod.names[i] <- paste0("(", i, ")")
}

tab1 <-
  stargazer(main_models_study1[[1]], main_models_study1[[2]], main_models_study1[[3]],
            omit = c("Age", "Female", "Education", "Income"), 
            type = "latex",
            title = "Study 1 Linear Models", 
            column.labels = mod.names,
            # type = "text",
            dep.var.labels = "Belief",
            add.lines=list(
              c("Demographic Controls", 
                rep("No", length(m1_s1)), 
                rep("No", length(m2_s1)), 
                rep("Yes", length(m3_s1))
              )
            ), 
            table.layout = "=ldc-ta-s-n",
            se = list(cl.rob.se.1[[1]], cl.rob.se.1[[2]], cl.rob.se.1[[3]],
                      cl.rob.se.2[[1]], cl.rob.se.2[[2]], cl.rob.se.2[[3]], cl.rob.se.2[[4]],
                      cl.rob.se.3[[1]], cl.rob.se.3[[2]], cl.rob.se.3[[3]], cl.rob.se.3[[4]]),
            p = list(p.vals1[[1]], p.vals1[[2]], p.vals1[[3]],
                     p.vals2[[1]], p.vals2[[2]], p.vals2[[3]], p.vals2[[4]],
                     p.vals3[[1]], p.vals3[[2]], p.vals3[[3]], p.vals3[[4]]),
            omit.stat = c("f", "ser"),
            # no.space = T,
            column.sep.width = "-5pt",
            intercept.bottom = FALSE,
            font.size = "tiny",
            order = c(1, 5),
            notes = c("Standard errors are clustered two-ways on Respondent and Narrative. ∗p<0.1;∗∗p<0.05;∗∗∗p<0.01"), 
            notes.align = "l",
            notes.append = FALSE,
            out.header = FALSE,
            header = FALSE
  ) %>% 
  star_insert_row(paste0("F Statistic", " ", 
                         "(df = ", 
                         #abs(cl.wald1[[1]][2,]$Df),  ";", " ", 
                         cl.wald1[[1]][2,]$Res.Df, ")", "&", 
                         round(cl.wald1[[1]][2,]$F, 3), "***", "&", 
                         round(cl.wald1[[2]][2,]$F, 3), "***", "&", 
                         round(cl.wald1[[3]][2,]$F, 3), "***", "&",
                         round(cl.wald2[[1]][2,]$F, 3), "***", "&",
                         round(cl.wald2[[2]][2,]$F, 3), "***", "&",
                         round(cl.wald2[[3]][2,]$F, 3), "***", "&",
                         round(cl.wald2[[4]][2,]$F, 3), "***", '&',
                         round(cl.wald3[[1]][2,]$F, 3), "***", "&",
                         round(cl.wald3[[2]][2,]$F, 3), "***", "&",
                         round(cl.wald3[[3]][2,]$F, 3), "***", "&",
                         round(cl.wald3[[4]][2,]$F, 3), "***",
                         "\\\\"), 
                  insert.after = 78) #%>% cat

write(tab1, file = "tables/tab1.tex")

# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Replicating Numeric and Non-Numeric CRT and 7 item vs 24 item AOT ####
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
m3_crt_s1 <- list()
m3_crt_s1[[1]] <- lm(response ~ CRT_numer*StoryTrue*RussianOrientation + 
                       Age + Female + factor(Education) + Income, psych_long)
m3_crt_s1[[2]] <- lm(response ~ CRT_nonnumer*StoryTrue*RussianOrientation + 
                       Age + Female + factor(Education) + Income, psych_long)
m3_crt_s1[[3]] <- lm(response ~ CRT_numer*StoryTrue*ProRussiaPoliOrientation + 
                       Age + Female + factor(Education) + Income, psych_long)
m3_crt_s1[[4]] <- lm(response ~ CRT_nonnumer*StoryTrue*ProRussiaPoliOrientation + 
                       Age + Female + factor(Education) + Income, psych_long)
m3_crt_s1[[5]] <- lm(response ~ AOT_short*StoryTrue*RussianOrientation + 
                       Age + Female + factor(Education) + Income, psych_long)
m3_crt_s1[[6]] <- lm(response ~ AOT_long*StoryTrue*RussianOrientation + 
                       Age + Female + factor(Education) + Income, psych_long)
m3_crt_s1[[7]] <- lm(response ~ AOT_short*StoryTrue*ProRussiaPoliOrientation + 
                       Age + Female + factor(Education) + Income, psych_long)
m3_crt_s1[[8]] <- lm(response ~ AOT_long*StoryTrue*ProRussiaPoliOrientation + 
                       Age + Female + factor(Education) + Income, psych_long)

vcov3 <- list()
vcov3 <- lapply(m3_crt_s1, cluster.vcov, ~ Respondent_Serial + question)

# get SEs
cl.rob.se.3 <- lapply(vcov3, function(x) sqrt(diag(x)))

# get F stats
cl.wald3 <- list()
for(i in 1:length(m3_crt_s1)) {
  cl.wald3[[i]]  <- lmtest::waldtest(m3_crt_s1[[i]], vcov = vcov3[[i]])
}

# get p-values
p.vals3 <- list()
for(i in 1:length(m3_crt_s1)) {
  p.vals3[[i]] <- lmtest::coeftest(m3_crt_s1[[i]], vcov = vcov3[[i]])[,4]
}

# exporting table reuslts
mod.names <- vector("character", length(m3_crt_s1))
for(i in 1:length(m3_crt_s1)) {
  mod.names[i] <- paste0("(", i, ")")
}

tab3 <- 
  stargazer(m3_crt_s1,
            omit = c("Age", "Female", "Education", "Income"),
            column.labels = mod.names,
            type = "latex",
            title = "Study 1 CRT Numeric vs Non-Numeric & AOT Short vs Long",
            # type = "text",
            dep.var.labels = "Belief",
            add.lines=list(c("Demographic Controls", rep("Yes", length(m3_crt_s1)))), 
            table.layout = "=ldc-ta-s-n",
            se = list( cl.rob.se.3[[1]], cl.rob.se.3[[2]], cl.rob.se.3[[3]], cl.rob.se.3[[4]],
                       cl.rob.se.3[[5]], cl.rob.se.3[[6]], cl.rob.se.3[[7]], cl.rob.se.3[[8]]),
            p = list( p.vals3[[1]], p.vals3[[2]], p.vals3[[3]], p.vals3[[4]],
                      p.vals3[[5]], p.vals3[[6]], p.vals3[[7]], p.vals3[[8]]),
            omit.stat = c("f", "ser"),
            # no.space = T,
            column.sep.width = "-5pt",
            intercept.bottom = FALSE,
            font.size = "tiny",
            order = c(1, 6),
            notes = c("Standard errors are clustered two-ways on Respondent and Narrative. ∗p<0.1;∗∗p<0.05;∗∗∗p<0.01"), 
            notes.align = "l",
            notes.append = FALSE,
            out.header = FALSE,
            header = FALSE
  ) %>% 
  star_insert_row(paste0("F Statistic", " ", 
                         "(df = ", 
                         #abs(cl.wald1[[1]][2,]$Df),  ";", " ", 
                         cl.wald3[[1]][2,]$Res.Df, ")", "&", 
                         round(cl.wald3[[1]][2,]$F, 3), "***", "&",
                         round(cl.wald3[[2]][2,]$F, 3), "***", "&",
                         round(cl.wald3[[3]][2,]$F, 3), "***", "&",
                         round(cl.wald3[[4]][2,]$F, 3), "***", '&',
                         round(cl.wald3[[5]][2,]$F, 3), "***", "&",
                         round(cl.wald3[[6]][2,]$F, 3), "***", "&",
                         round(cl.wald3[[7]][2,]$F, 3), "***", "&",
                         round(cl.wald3[[8]][2,]$F, 3), "***",
                         "\\\\"), 
                  insert.after = 108) # %>% cat

write(tab3, file = "tables/tab3.tex")

# ++++++++++++++++++++++++++++++++++++++++++++++
# Appendix Theme and Strategy Fixed Effects ####
# ++++++++++++++++++++++++++++++++++++++++++++++

# need to create dummies so that we have anti-west baseline
psych_long$StrategyUndermineUkraine <- ifelse(psych_long$strategy == 3, 1, 0)
psych_long$StrategyProRussia <- ifelse(psych_long$strategy == 2, 1, 0)

# part 4 (with FEs for theme + strategy)
m4_s1a <- list()
m4_s1a[[1]] <- lm(response ~ CRT*StoryTrue*RussianOrientation + Age + Female + 
                    factor(Education) + Income + as_factor(theme) + 
                    StrategyUndermineUkraine + StrategyProRussia, psych_long)
m4_s1a[[2]] <- lm(response ~ AOT*StoryTrue*RussianOrientation + Age + Female + 
                    factor(Education) + Income + as_factor(theme) + 
                    StrategyUndermineUkraine + StrategyProRussia, psych_long)
m4_s1a[[3]] <- lm(response ~ CRT*StoryTrue*ProRussiaPoliOrientation + Age + Female + 
                    factor(Education) + Income + as_factor(theme) + 
                    StrategyUndermineUkraine + StrategyProRussia, psych_long)
m4_s1a[[4]] <- lm(response ~ AOT*StoryTrue*ProRussiaPoliOrientation + Age + Female + 
                    factor(Education) + Income + as_factor(theme) + 
                    StrategyUndermineUkraine + StrategyProRussia, psych_long)

vcov4a <- list()
vcov4a <- lapply(m4_s1a, cluster.vcov, ~ Respondent_Serial + question)

# get SEs
cl.rob.se.4a <- lapply(vcov4a, function(x) sqrt(diag(x)))

# get F-stats
cl.wald4a <- list()
for(i in 1:length(m4_s1a)) {
  cl.wald4a[[i]]  <- lmtest::waldtest(m4_s1a[[i]], vcov = vcov4a[[i]])
}

# get p-values
p.vals4 <- list()
for(i in 1:length(m4_s1a)) {
  p.vals4[[i]] <- lmtest::coeftest(m4_s1a[[i]], vcov = vcov4a[[i]])[,4]
}

# exporting table results
mod.names <- vector("character", length(m4_s1a))
for(i in 1:length(m4_s1a)) {
  mod.names[i] <- paste0("(", i, ")")
}

reg_m4_s1a <-  # strategy baseline = anti-west
  stargazer(m4_s1a,
            #omit = c("Age", "Female", "Education", "Income", "StrategyUndermineUkraine", "StrategyProRussia", "theme"), 
            type = "latex",
            title = "Study 1 Linear Models with Theme and Strategy Controls", 
            column.labels = mod.names,
            # type = "text",
            dep.var.labels = "Belief",
            # add.lines=list(c("Demographic Controls", rep("Yes", length(m4_s1a))),
            #                c("Theme and Strategy Controls", rep("Yes", length(m4_s1a)))), 
            table.layout = "=ldc-ta-s-n",
            se = list( cl.rob.se.4a[[1]], cl.rob.se.4a[[2]], cl.rob.se.4a[[3]], cl.rob.se.4a[[4]]),
            p = list( p.vals4[[1]], p.vals4[[2]], p.vals4[[3]], p.vals4[[4]]),
            omit.stat = c("f", "ser"),
            # no.space = T,
            column.sep.width = "-5pt",
            intercept.bottom = FALSE,
            font.size = "tiny",
            order = c(1, 6),
            notes = "This will be replaced",
            notes.align = "l",
            notes.append = FALSE,
            out.header = FALSE,
            header = FALSE
  ) %>% 
  star_insert_row(paste0("F Statistic", " ", 
                         "(df = ", 
                         #abs(cl.wald1[[1]][2,]$Df),  ";", " ", 
                         cl.wald4a[[1]][2,]$Res.Df, ")", "&", 
                         round(cl.wald4a[[1]][2,]$F, 3), "***", "&",
                         round(cl.wald4a[[2]][2,]$F, 3), "***", "&",
                         round(cl.wald4a[[3]][2,]$F, 3), "***", "&",
                         round(cl.wald4a[[4]][2,]$F, 3), "***",
                         "\\\\"), 
                  # insert.after = 73) # with omitting coefficients
                  insert.after = 110) # no omitting coefficients

note.latex <- "\\multicolumn{5}{l} {\\parbox[t]{11cm}{ \\textit{Notes:} Standard errors are clustered two-ways on Respondent and Narrative. ∗p<0.1;∗∗p<0.05;∗∗∗p<0.0}} \\\\"
reg_m4_s1a[grepl("Note",reg_m4_s1a)] <- note.latex
cat (reg_m4_s1a, sep = "\n")

write(reg_m4_s1a, file = "tables/reg4_s1a.tex")

rm(psych_long)

# +++++++++++++++++++++++++++++++++++++++++++++++++
# Un-Standardized (unscaled) Measures Analyses ####
# +++++++++++++++++++++++++++++++++++++++++++++++++

psych_long_us <- read_dta("./data_clean/psych_clean_long.dta") %>% 
  rename(StoryTrue = true,
         # CRT = zCRT_all,
         CRT = CRT_all, #CRT IS NOW UNSCALED VERSION
         # AOT = AOT_all,
         AOT =  AOT_allUnscaled,
         # AOT_shortUnscaled,
         # AOT_longUnscaled,
         # RussianOrientation = russianOrientation,
         # ProRussiaPoliOrientation = zantiEuropeOrientation,
         RussianOrientation = XrussianOrientation,
         ProRussiaPoliOrientation = antiEuropeOrientation,
         Female = RSPSEX,
         Age = RSPAGE,
         Education = RSPEDUC,
         Income = income) %>% 
  mutate(Belief = response,
         AOT = (AOT - 1)/4)

# Linear Models

# part 1
m1_s1_us <- list()
m1_s1_us[[1]] <- lm(response ~ CRT*StoryTrue, psych_long_us)
m1_s1_us[[2]] <- lm(response ~ AOT*StoryTrue, psych_long_us)
m1_s1_us[[3]] <- lm(response ~ cogStyleAgg*StoryTrue, psych_long_us)

vcov1_s1_us <- list()
cl <- makeCluster(4)
vcov1_s1_us <- lapply(m1_s1_us, cluster.vcov, ~ Respondent_Serial + question, parallel = TRUE)
stopCluster(cl)

# part 2
m2_s1_us <- list()
m2_s1_us[[1]] <- lm(response ~ CRT*StoryTrue*RussianOrientation, psych_long_us)
m2_s1_us[[2]] <- lm(response ~ AOT*StoryTrue*RussianOrientation, psych_long_us)
m2_s1_us[[3]] <- lm(response ~ CRT*StoryTrue*ProRussiaPoliOrientation, psych_long_us)
m2_s1_us[[4]] <- lm(response ~ AOT*StoryTrue*ProRussiaPoliOrientation, psych_long_us)

cl <- vcov2_s1_us <- list()
cl <- makeCluster(4)
vcov2_s1_us <- lapply(m2_s1_us, cluster.vcov, ~ Respondent_Serial + question, parallel = TRUE)
stopCluster(cl)

# saving main models as .rds
write_rds(m2_s1_us, file = "./data_clean/study1_main_models_unscaled.rds")


# part 3
m3_s1_us <- list()
m3_s1_us[[1]] <- lm(response ~ CRT*StoryTrue*RussianOrientation + Age + Female + factor(Education) + Income, 
                    psych_long_us)
m3_s1_us[[2]] <- lm(response ~ AOT*StoryTrue*RussianOrientation + Age + Female + factor(Education) + Income, 
                    psych_long_us)
m3_s1_us[[3]] <- lm(response ~ CRT*StoryTrue*ProRussiaPoliOrientation + Age + Female + factor(Education) + Income, 
                    psych_long_us)
m3_s1_us[[4]] <- lm(response ~ AOT*StoryTrue*ProRussiaPoliOrientation + Age + Female + factor(Education) + Income, 
                    psych_long_us)

vcov3_s1_us <- list()
cl <- makeCluster(4)
vcov3_s1_us <- lapply(m3_s1_us, cluster.vcov, ~ Respondent_Serial + question, parallel = TRUE)
stopCluster(cl)

# part 4 (with FEs for theme)
m4_s1_us <- list()
m4_s1_us[[1]] <- lm(response ~ CRT*StoryTrue*RussianOrientation + Age + Female + factor(Education) + Income +
                      factor(theme), 
                    psych_long_us)
m4_s1_us[[2]] <- lm(response ~ AOT*StoryTrue*RussianOrientation + Age + Female + factor(Education) + Income + 
                      factor(theme), 
                    psych_long_us)
m4_s1_us[[3]] <- lm(response ~ CRT*StoryTrue*ProRussiaPoliOrientation + Age + Female + factor(Education) + Income + 
                      factor(theme), 
                    psych_long_us)
m4_s1_us[[4]] <- lm(response ~ AOT*StoryTrue*ProRussiaPoliOrientation + Age + Female + factor(Education) + Income + 
                      factor(theme), 
                    psych_long_us)

vcov4_s1_us <- list()
cl <- makeCluster(4)
vcov4_s1_us <- lapply(m4_s1_us, cluster.vcov, ~ Respondent_Serial + question, parallel = TRUE)
stopCluster(cl)

# get SEs
cl.rob.se.1 <- lapply(vcov1_s1_us, function(x) sqrt(diag(x)))
cl.rob.se.2 <- lapply(vcov2_s1_us, function(x) sqrt(diag(x)))
cl.rob.se.3 <- lapply(vcov3_s1_us, function(x) sqrt(diag(x)))

# get F Stats
cl.wald1 <- list()
for(i in 1:length(m1_s1_us)) {
  cl.wald1[[i]]  <- waldtest(m1_s1_us[[i]], vcov = vcov1_s1_us[[i]])
}

cl.wald2 <- list()
for(i in 1:length(m2_s1_us)) {
  cl.wald2[[i]]  <- waldtest(m2_s1_us[[i]], vcov = vcov2_s1_us[[i]])
}

cl.wald3 <- list()
for(i in 1:length(m3_s1_us)) {
  cl.wald3[[i]]  <- waldtest(m3_s1_us[[i]], vcov = vcov3[[i]])
}

# get p-values
p.vals1 <- list()
for(i in 1:length(m1_s1_us)) {
  p.vals1[[i]] <- coeftest(m1_s1_us[[i]], vcov = vcov1_s1_us[[i]])[,4]
}

p.vals2 <- list()
for(i in 1:length(m2_s1_us)) {
  p.vals2[[i]] <- coeftest(m2_s1_us[[i]], vcov = vcov2_s1_us[[i]])[,4]
}

p.vals3 <- list()
for(i in 1:length(m3_s1_us)) {
  p.vals3[[i]] <- coeftest(m3_s1_us[[i]], vcov = vcov3_s1_us[[i]])[,4]
}

# saving main models
main_models_us <- list(m1_s1_us, m2_s1_us, m3_s1_us)
save(main_models_us, file =  "./data_clean/main_clustered_models_study1_unscaled.Rdata")

# exporting table results
mod.names <- vector("character", sum(length(m1_s1_us), length(m2_s1_us), length(m3_s1_us)))
for(i in 1:sum(length(m1_s1_us), length(m2_s1_us), length(m3_s1_us))) {
  mod.names[i] <- paste0("(", i, ")")
}

tab1 <-
  stargazer(m1_s1_us, m2_s1_us, m3_s1_us,
            omit = c("Age", "Female", "Education", "Income"), 
            type = "latex",
            title = "Study 1 Linear Models", 
            column.labels = mod.names,
            # type = "text",
            dep.var.labels = "Belief",
            add.lines=list(
              c("Demographic Controls", 
                rep("No", length(m1_s1_us)), 
                rep("No", length(m2_s1_us)), 
                rep("Yes", length(m3_s1_us))
              )
            ), 
            table.layout = "=ldc-ta-s-n",
            se = list(cl.rob.se.1[[1]], cl.rob.se.1[[2]], cl.rob.se.1[[3]],
                      cl.rob.se.2[[1]], cl.rob.se.2[[2]], cl.rob.se.2[[3]], cl.rob.se.2[[4]],
                      cl.rob.se.3[[1]], cl.rob.se.3[[2]], cl.rob.se.3[[3]], cl.rob.se.3[[4]]),
            p = list(p.vals1[[1]], p.vals1[[2]], p.vals1[[3]],
                     p.vals2[[1]], p.vals2[[2]], p.vals2[[3]], p.vals2[[4]],
                     p.vals3[[1]], p.vals3[[2]], p.vals3[[3]], p.vals3[[4]]),
            omit.stat = c("f", "ser"),
            # no.space = T,
            column.sep.width = "-5pt",
            intercept.bottom = FALSE,
            font.size = "tiny",
            order = c(1, 5),
            notes = c("Standard errors are clustered two-ways on Respondent and Narrative. ∗p<0.1;∗∗p<0.05;∗∗∗p<0.01"), 
            notes.align = "l",
            notes.append = FALSE,
            out.header = FALSE,
            header = FALSE
  ) %>% 
  star_insert_row(paste0("F Statistic", " ", 
                         "(df = ", 
                         #abs(cl.wald1[[1]][2,]$Df),  ";", " ", 
                         cl.wald1[[1]][2,]$Res.Df, ")", "&", 
                         round(cl.wald1[[1]][2,]$F, 3), "***", "&", 
                         round(cl.wald1[[2]][2,]$F, 3), "***", "&", 
                         round(cl.wald1[[3]][2,]$F, 3), "***", "&",
                         round(cl.wald2[[1]][2,]$F, 3), "***", "&",
                         round(cl.wald2[[2]][2,]$F, 3), "***", "&",
                         round(cl.wald2[[3]][2,]$F, 3), "***", "&",
                         round(cl.wald2[[4]][2,]$F, 3), "***", '&',
                         round(cl.wald3[[1]][2,]$F, 3), "***", "&",
                         round(cl.wald3[[2]][2,]$F, 3), "***", "&",
                         round(cl.wald3[[3]][2,]$F, 3), "***", "&",
                         round(cl.wald3[[4]][2,]$F, 3), "***",
                         "\\\\"), 
                  insert.after = 78) #%>% cat

write(tab1, file = "tables/tab1_unscaled.tex")

# HTML version
tab2 <-
  stargazer(m1_s1_us, m2_s1_us, m3_s1_us,
            omit = c("Age", "Female", "Education", "Income"), 
            type = "html",
            title = "Study 1 Linear Models", 
            column.labels = mod.names,
            # type = "text",
            dep.var.labels = "Belief",
            add.lines=list(
              c("Demographic Controls", 
                rep("No", length(m1_s1_us)), 
                rep("No", length(m2_s1_us)), 
                rep("Yes", length(m3_s1_us))
              )
            ), 
            table.layout = "=ldc-ta-s-n",
            se = list(cl.rob.se.1[[1]], cl.rob.se.1[[2]], cl.rob.se.1[[3]],
                      cl.rob.se.2[[1]], cl.rob.se.2[[2]], cl.rob.se.2[[3]], cl.rob.se.2[[4]],
                      cl.rob.se.3[[1]], cl.rob.se.3[[2]], cl.rob.se.3[[3]], cl.rob.se.3[[4]]),
            p = list(p.vals1[[1]], p.vals1[[2]], p.vals1[[3]],
                     p.vals2[[1]], p.vals2[[2]], p.vals2[[3]], p.vals2[[4]],
                     p.vals3[[1]], p.vals3[[2]], p.vals3[[3]], p.vals3[[4]]),
            omit.stat = c("f", "ser"),
            # no.space = T,
            column.sep.width = "-5pt",
            intercept.bottom = FALSE,
            font.size = "tiny",
            order = c(1, 5),
            notes = c("Standard errors are clustered two-ways on Respondent and Narrative. ∗p<0.1;∗∗p<0.05;∗∗∗p<0.01"), 
            notes.align = "l",
            notes.append = FALSE,
            out.header = FALSE,
            header = FALSE
  ) %>% 
  star_insert_row(paste0(
    "<tr><td style=\"text-align:left\">",
    "F Statistic", " ", "(df = ", cl.wald1[[1]][2,]$Res.Df, ")", "</td><td>", 
    round(cl.wald1[[1]][2,]$F, 3), "***", "</td><td>", 
    round(cl.wald1[[2]][2,]$F, 3), "***", "</td><td>", 
    round(cl.wald1[[3]][2,]$F, 3), "***", "</td><td>",
    round(cl.wald2[[1]][2,]$F, 3), "***", "</td><td>",
    round(cl.wald2[[2]][2,]$F, 3), "***", "</td><td>",
    round(cl.wald2[[3]][2,]$F, 3), "***", "</td><td>",
    round(cl.wald2[[4]][2,]$F, 3), "***", '</td><td>',
    round(cl.wald3[[1]][2,]$F, 3), "***", "</td><td>",
    round(cl.wald3[[2]][2,]$F, 3), "***", "</td><td>",
    round(cl.wald3[[3]][2,]$F, 3), "***", "</td><td>",
    round(cl.wald3[[4]][2,]$F, 3), "***", "</td></tr>"#,
    # "</table>"
    #"\\\\"
    ), 
    insert.after = 78) #%>% cat

# need to manually go into html to insert F-stat in right spot

cm <- c('(Intercept)' = 'Constant',
        'StoryTrue' = 'True',
        'CRT' = 'CRT',
        'AOT' = 'AOT',
        'RussianOrientation' = 'Russian Orientation',
        'ProRussiaPoliOrientation' = 'Russia Pol. Orientation',
        'CRT:StoryTrue' = 'CRT×True',
        'AOT:StoryTrue' = 'AOT×True',
        'CRT:RussianOrientation'= 'CRT×RussianOrientation',
        'AOT:RussianOrientation' = 'AOT×RussianOrientation',
        'StoryTrue:RussianOrientation' = 'True×RussianOrientation',
        'CRT:StoryTrue:RussianOrientation'    = 'CRT×True×RussianOrientation',
        'AOT:StoryTrue:RussianOrientation'    = 'AOT×True×RussianOrientation',
        'CRT:ProRussiaPoliOrientation' = 'CRT×Pro Russia Pol. Orientation',
        'AOT:ProRussiaPoliOrientation'  = 'AOT×Pro Russia Pol. Orientation',
        'StoryTrue:ProRussiaPoliOrientation' = 'True×Pro Russia Pol. Orientation',
        'CRT:StoryTrue:ProRussiaPoliOrientation' = 'CRT×True×Pro Russia Pol. Orientation',
        'AOT:StoryTrue:ProRussiaPoliOrientation' = 'AOT×True×Pro Russia Pol. Orientation')


f <- function(x) format(round(x, 3), big.mark=",")
gm <- list(
  list("raw" = "nobs", "clean" = "Observations", "fmt" = f),
  list("raw" = "adj.r.squared", "clean" = "Adjusted R^2", "fmt" = f))

msummary(c(m2_s1_us), output = "./tables/study_1_main_paper.html",
         title = "Study 1: Linear Models",
         notes	= c("Two-way standard errors clustered on Respondent and Narrative"),
         statistic_override = c(vcov2_s1_us),
         coef_map = cm,
         gof_map = gm,
         escape = TRUE,
         stars = TRUE)   #tab_spanner(label = 'Literacy', columns = c('OLS 1', 'NBin 1')) %>% # add columns


write(tab2, file = "./tables/tab1_unscaled.html")

